/*
 * MineralKinetics.cpp
 *
 *  Created on: 21 Jul 2011
 *      Author: Allan
 */

#include "MineralKinetics.h"

// C++ includes
#include <cmath>
using namespace std;

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "Core/Multiphase.h"
#include "Utils/Assert.h"

MineralKineticsInfo::MineralKineticsInfo() 
{}

MineralKineticsInfo::MineralKineticsInfo(const string& mineral) :
mineral(mineral)
{}

void MineralKineticsInfo::AddMineralMechanism(double k, double Ea) 
{
	MineralMechanism mechanism;
	
	mechanism.k  = k;
	mechanism.Ea = Ea;
	
	mechanisms.push_back(mechanism);
}

void MineralKineticsInfo::AddMineralMechanism(double k, double Ea, const vector<tuple<string,string,double>>& catalysts) 
{
	MineralMechanism mechanism;
	
	mechanism.k  = k;
	mechanism.Ea = Ea;
	
	BOOST_FOREACH(auto c, catalysts)
	{
		MineralCatalyst catalyst;
		
		catalyst.name     = get<0>(c);
		catalyst.quantity = get<1>(c);
		catalyst.power    = get<2>(c);
		
		mechanism.catalysts.push_back(catalyst);
	}
	
	mechanisms.push_back(mechanism);
}

void MineralKineticsInfo::SetSpecificSurfaceArea(double ssa)
{
	specificSurfaceArea = ssa;
}

const string MineralKineticsInfo::GetMineralName() const 
{
	return mineral;
}

const double MineralKineticsInfo::GetSpecificSurfaceArea() const 
{
	return specificSurfaceArea; 
}

const vector<MineralMechanism> MineralKineticsInfo::GetMineralMechanisms() const 
{
	return mechanisms;
}

MineralKinetics::MineralKinetics() :
imineral(-1), specificSurfaceArea(0.0)
{}

MineralKinetics::MineralKinetics(const MineralKineticsInfo& info, const Multiphase& multiphase) :
mineral(info.GetMineralName()), mechanisms(info.GetMineralMechanisms()), 
imineral(multiphase[mineral]), specificSurfaceArea(info.GetSpecificSurfaceArea())
{
	// Set the indexes of the catalysts species in the mineral mechanisms
	BOOST_FOREACH(MineralMechanism& mechanism, mechanisms)
		BOOST_FOREACH(MineralCatalyst& catalyst, mechanism.catalysts)
			catalyst.index = multiphase[catalyst.name];
}

MineralKinetics::~MineralKinetics()
{}

//const double MineralKinetics::KineticRate(const Reaction& reaction, double T, double P, const VectorXd& n, const VectorXd& a) const
//{
//	// The universal gas constant in kJ/(mol.K)
//	const double R = 8.314472E-3;
//	
//	// The kinetic rate constant for the neutral mechanism
//	const double k = mechanisms[0].k;
//	
//	// The Arrhenius activation energy for the neutral mechanism
//	const double Ea = mechanisms[0].Ea;
//	
//	// The reactive surface area of the mineral, which is given by the product of the mineral specific surface area and its current number of moles
//	const double A = specificSurfaceArea * n[imineral];
//	
//	// The equilibrium constant of the mineral reaction
//	const double K = reaction.EquilibriumConst(T, P);
//	
//	// The reaction quotient of the mineral reaction
//	const double Q = reaction.ReactionQuotient(a);
//	
//	// Calculate the rate of dissolution/precipitation of the mineral (positive: dissolution, negative: precipitation)
//	return k * exp(-Ea/R * (1.0/(T + 273.15) - 1.0/298.15)) * A * (1.0 - Q/K);
//}

const double MineralKinetics::KineticRate(const Reaction& reaction, double T, double P, const VectorXd& n, const VectorXd& a) const
{
	// The temperature in Kelvin
	const double Tk = T + 273.15;
	
	// The universal gas constant in kJ/(mol.K)
	const double R = 8.314472E-3;
	
	// The equilibrium constant of the mineral reaction
	const double K = reaction.EquilibriumConst(T, P);
	
	// The reaction quotient of the mineral reaction
	const double Q = reaction.ReactionQuotient(a);
	
	// The reactive surface area of the mineral, which is given by the product of the mineral specific surface area and its current number of moles
	const double A = specificSurfaceArea * n[imineral];
	
	// The effective rate constant taking into account all mechanisms
	double keff = 0.0;
	
	// Iterate over all the mechanisms
	BOOST_FOREACH(const MineralMechanism& mechanism, mechanisms)
	{
		// The kinetic rate constant for the current mechanism at 25 C
		const double k25 = mechanism.k;
		
		// The Arrhenius activation energy for the current mechanism
		const double Ea = mechanism.Ea;
		
		// Calculate the rate constant for the current mechanism
		double k = k25 * exp(-Ea/R * (1.0/Tk - 1.0/298.15));
		
		// Update the rate constant for the current mechanism with the catalyst contributions
		BOOST_FOREACH(const MineralCatalyst& catalyst, mechanism.catalysts)
			k *= pow(a[catalyst.index], catalyst.power);
		
		// Update the effective rate constant taking into account all mechanisms
		keff += k;
	}
	
	// Calculate the rate of dissolution/precipitation of the mineral (positive: dissolution, negative: precipitation)
	return keff * A * (1.0 - Q/K);
}

const VectorXd MineralKinetics::DerivativeKineticRate(const Reaction& reaction, double T, double P, const VectorXd& n, const VectorXd& a) const
{
	return VectorXd::Zero(n.rows());
}
